program Euler_Heun
implicit none
  integer         :: i
  real            :: x1,x2,y1,y2
  real,parameter  :: h = .5

x1 = 0
y1 = 1
! implementing Euler's method
open(1,file='euler')
write(1,*)x1,y1

do i= 1,8

y2=y1+h*f(x1,y1)

x2 = x1+0.25
x1 = x2
y1 = y2
write(1,*)x2,y2
end do

!implementing Heun's method
x1 = 0
y1 = 1

open(2,file='heun')
write(2,*)x1,y1

do i= 1,8

y2=y1+h*(f(x1,y1)+f(x1+h,y1+h*f(x1,y1)))/2

x2 = x1+0.25
x1 = x2
y1 = y2
write(2,*)x2,y2
end do


contains
real function f(a,b)
real,intent(in) :: a,b
f = (1+2*a)*sqrt(b)
end function
end program 

